all=importdata('D:\study\bsegit\bse-pap\hanchen_first_neutron1124.out');
tspan= importdata('all_t.xls' );

% all=load('D:\study\bsegit\MW\hmxb.dat');
% 1  2   3  4   5  6   7   8   9  10 11 12  13    14 15 16  17  18 19
% m1,m2,Lx,ecc0,a,tb0,mt1,mt2,kw,kw2,mx,mx2,eccx,tbx,f,nnn,uuu,m_tr,i
% m1,m2,Lx,ecc0,a,tb0,mt1,mt2,kw,kw2,mx,mx2,eccx,tbx,f,m_transferrate,b_iso,m_tr,i
% 20 \delta T
global  AllmassNUM M1MAX M1MIN SFR qNUM aNUM alpha

alpha=2.7;
tmax=1d4*1d6; % in units of year
AllmassNUM=1d3;
M1MAX=150;
M1MIN=2;
SFR=GetSFR(43.5,alpha,0.535,0,M1MIN,5);
% SFR=10.6984;
qNUM=9;
aNUM=100;
yearsc=3.1557d+07;
Msun=1.9891d30;
yearday=365.0;

all_1=all;
m_loc=4;
tb_loc=6;

m_loc=11;
tb_loc=14;

mstep=0.1;
% mline=[9.625:mstep:20.625,21.25:0.5:30.25];
mline = 2-0.125 : 0.25  : 29.5+ 0.125;
tbline=(0.2 - 0.05 : 0.1 : 3.0 + 0.05);
% mline=[0.1:mstep:20.625];
% tbline=(0.001:0.1:3.0);

% change orbit  to circle
ecc = all_1(:, 13);
all_1(:,tb_loc) = all_1(:,tb_loc) .*yearday .* (1 - ecc.^2 ).^ 1.5;

all_1(all_1(:,m_loc)<=min(mline),:)=[];
all_1(all_1(:,m_loc)>= max(mline),:)=[];
all_1(all_1(:,tb_loc)<= 10.^(min(tbline)),:)=[];
all_1(all_1(:,tb_loc)>= 10.^(max(tbline)),:)=[];
% all_1(all_1(:,9)>2,:)=[];
% all_1(all_1(:,15)>100,:)=[];

%check shaoli
plot(all_1(:,m_loc),log10(all_1(:,tb_loc)),'.')


[counts,centers] = hist(log10(all_1(:,tb_loc)));
bar(centers,log10(counts))
title('')
xlabel('log10(P_{orb} (day))')
ylabel('log10(N)')
ylim([2 4])
[counts,centers] = hist((all_1(:,m_loc)));
bar(centers,log10(counts))
title('')
xlabel('log10(P_{orb} (day))')
ylabel('log10(N)')
ylim([2 4])
mpos=sum(all_1(:,m_loc)>mline,2);
tbpos=sum(all_1(:,tb_loc)>10.^tbline,2);
[len,~]=size(all_1);
% result=zeros(length(mline),length(tbline));
result=zeros(length(tbline),length(mline));
for i=1:len
    result(tbpos(i),mpos(i))=result(tbpos(i),mpos(i))+...
        massFunction(all_1(i,1));
%     result(tbpos(i),mpos(i))=result(tbpos(i),mpos(i))+...
%         1;
end

a=pcolor(mline,tbline,result);
a.EdgeColor='none';

result = result(1:end-1, 1:end-1 );


